The Mitchell River is a river located in Far North Queensland, Australia. The river rises on the Atherton Tableland about 50 kilometres (31 mi) northwest of Cairns, and flows about 750 kilometres (470 mi) northwest across Cape York Peninsula from Mareeba to the Gulf of Carpentaria. we will use this case study to present some functionalities of the prioriactions package.
We started loading libraries:
library(prioritizr) #To use the simulate features distribution
library(prioriactions)
library(raster) #To plot of shapefiles
library(tmap) #To create cool maps
library(scales) #To standardize the value of amount
library(reshape2) #To use the melt function
library(sp) #To use the spplot function
library(viridis) #To use viridis pallete
We divided the whole catchment (71,630 km2 into 2316 sites (i.e., sub-catchments), each one included the portion of river length between two consecutive river connections. We considered four major threats to freshwater fish species in the catchment: water buffalo (Bubalis bubalis), cane toad (Bufo marinus), river flow alteration (caused by impoundments, channels for water extractions and levee banks) and grazing land use. Also, we used the modelled spatial distribution of 45 fish species in the Mitchell river catchment as our conservation features. The distribution of features and threats () will be simulated while the other input files a prioriactions will be loaded directly from the package as follows:
path <- system.file("extdata/input_big/", package = "prioriactions")
pu_data <- data.table::fread(paste0(path,"unitCost_Big.csv"), data.table = FALSE)
features_data <- data.table::fread(paste0(path,"target_Big.csv"), data.table = FALSE)
bound_data <- data.table::fread(paste0(path,"boundary_Big.csv"), data.table = FALSE)
sensitivity_data <- data.table::fread(paste0(path,"sensibility_Big.csv"), data.table = FALSE)
We load the shapefile of the case study also included as part of the package installation:
#reading shapefile
shp_mitchell = raster::shapefile(paste0(path,"Mitchell.shp"))
sp::spplot(shp_mitchell, zcol = "Shape_Area", names.attr = "Area",
main = "Planning unit areas", col.regions = viridis::viridis(20))
To model the distribution of features and threats, we use the function simulate_species from prioritizr package. Due this simulation is performed on a raster object, we rasterize the shapefile (through the rasterize function of the raster library) and then vectorize it again with the extract function of the same library.
#We set the seed to can replicate the same simulations
set.seed(1)
#We transform the shapefile to raster with a convenient resolution for simulation
ext <- raster::extent(shp_mitchell)
r <- raster::raster(ext, res = 1000)
r <- raster::rasterize(shp_mitchell, r, field=1)
#Simulation of 45 features distribution
dist_features <- prioritizr::simulate_species(r, 45, model = RandomFields::RMgauss(scale = 40000))
#> .............................................
#We project the simulated rasters to shapefile
rij_data <- raster::extract(dist_features, shp_mitchell, fun = mean)
colnames(rij_data) <- features_data$id
In this way we obtain spatially auto correlated features distributions. Thus, for example, the distribution of simulated feature 1 is:
shp_mitchell$specie_distribution <- rij_data[,1]
#Plot distribution with tmap library
tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("specie_distribution", pal = c("white", "seagreen")) +
tmap::tm_borders(col="black", lwd = 0.5)
Instead, the distribution of simulated threat 1 is:
#Simulation of 4 threats distribution
dist_threats <- prioritizr::simulate_species(r, 4, model = RandomFields::RMgauss(scale = 40000))
#> ....
#We project the simulated rasters to shapefile
threats_data <- raster::extract(dist_threats, shp_mitchell, fun = mean)
colnames(threats_data) <- 1:4
shp_mitchell$threat_distribution <- threats_data[,1]
#Plot distribution with tmap library
tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("threat_distribution", pal = c("white", "red4")) +
tmap::tm_borders(col="black", lwd = 0.5)
For the first prioritization exercise, we transform the amount data to values 0 and 1 for both features distribution and threats. In the case of features, a threshold of 0.5 was set, that is, if the value of the amount is equal to or greater than 0.5, it is considered that the features is present at the site. Instead, a threat is considered present if it has any non-zero amount value.
#Modifying continuous to binary amount data of features
rij_data_binary <- ifelse(rij_data > 0.5, 1, 0)
shp_mitchell$feature_distribution <- rij_data_binary[,1]
#Modifying continuous to binary amount data of threats
threats_data_binary <- ifelse(threats_data > 0.0, 1, 0)
#Plot example of feature distribution with tmap library
tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("feature_distribution", pal = c("white", "seagreen"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
Then, we create the features and threat distribution inputs in the corresponding format that will be used by the prioriactions package. Also, we modify the features_data to set the target of each feature at 15%.
#Feature distribution data-------------------------------------
#We get the extended matrix of distribution data for both cases, set the specific col names and we are left with only positive values
rij_data_binary_large <- reshape2::melt(rij_data_binary)
colnames(rij_data_binary_large) <- c("pu", "species", "amount")
rij_data_binary_large <- rij_data_binary_large[rij_data_binary_large$amount > 0, ]
#threat distribution data-------------------------------------
threats_data_binary_large <- reshape2::melt(threats_data_binary)
colnames(threats_data_binary_large) <- c("pu", "threats", "amount")
threats_data_binary_large <- threats_data_binary_large[threats_data_binary_large$amount > 0, ]
threats_data_binary_large$cost <- 1
threats_data_binary_large$status <- 0
#features data-------------------------------------
features_data$target <- colSums(rij_data_binary)*0.15
The base model considers the prioritization of conservation actions using the presence / absence of features and threats. In addition, a target has been set for all features of 15% of their representatives in the basin. The model, for its part, does not consider spatial requirements to the actions (blm and blm_actions parameters equal to 0. Furthermore, there is no attempt to optimize that planning actions occur within the same planning units (i.e., the curve parameter is equal to 1). The solver was configured to stop when a gap of at least 5% is achieved (0% meaning that at least one of the optimal solutions has been found).
input_data.base <- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_binary_large,
threats = threats_data_binary_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Registered S3 methods overwritten by 'prioriactions':
#> method from
#> print.ConservationProblem prioritizr
#> print.OptimizationProblem prioritizr
#> Correctly loaded inputs
model.base <- prioriactions::min_costs(input_data.base, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.base <- prioriactions::solve(model.base, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 106792 nonzeros
#> Model fingerprint: 0x3d759f11
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#> Matrix range [3e-01, 4e+00]
#> Objective range [1e+00, 1e+00]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+02, 2e+02]
#> Found heuristic solution: objective 3304.0000000
#> Found heuristic solution: objective 2145.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.07s
#> Presolved: 2361 rows, 11545 columns, 106757 nonzeros
#> Variable types: 0 continuous, 11545 integer (11545 binary)
#>
#> Root relaxation: objective 8.769123e+02, 3052 iterations, 0.19 seconds
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 876.91232 0 694 2145.00000 876.91232 59.1% - 0s
#> H 0 0 1397.0000000 876.91232 37.2% - 0s
#> H 0 0 1396.0000000 876.91232 37.2% - 0s
#> 0 0 1029.85774 0 840 1396.00000 1029.85774 26.2% - 0s
#> H 0 0 1362.0000000 1029.85774 24.4% - 0s
#> 0 0 1036.03426 0 885 1362.00000 1036.03426 23.9% - 1s
#> 0 0 1134.80250 0 739 1362.00000 1134.80250 16.7% - 1s
#>
#> Cutting planes:
#> Cover: 1352
#> Clique: 10
#> MIR: 36
#>
#> Explored 1 nodes (7375 simplex iterations) in 1.38 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 5: 1362 1396 1397 ... 3304
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.362000000000e+03, best bound 1.135000000000e+03, gap 16.6667%
Note that we have achieved a gap of 0.24% and a objective value of 1265. To obtain the distribution of conservation actions we use the getSolutionActions() function, while that to obtain the vector of sites that are incorporated into the conservation plan we use the getSolutionUnits() function. Note that planning units may be selected just for adding connectivity to the plan without taking action within these sites.
#Getting unit distribution selected
solution_units.base <- solution.base$getSolutionUnits()
#Assign solution to shapefile field to plot it
shp_mitchell$sol.base <- solution_units.base$solution
#Plot of action 1 distribution
tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sol.base", pal = c("white", "gray40"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
Whereas, if we want to obtain the actions to be carried out within the sites:
#Getting action distribution
solution_actions.base <- solution.base$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.base <- solution_actions.base$`1`
shp_mitchell$action_2.base <- solution_actions.base$`2`
shp_mitchell$action_3.base <- solution_actions.base$`3`
shp_mitchell$action_4.base <- solution_actions.base$`4`
#Actions plots
plot_action1.base <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_1.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action2.base <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_2.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action3.base <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_3.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action4.base <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_4.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
tmap::tmap_arrange(plot_action1.base , plot_action2.base , plot_action3.base , plot_action4.base)
Or the distribution of the sum of the actions (higher density of actions):
shp_mitchell$sum_actions.base <- solution_units.base$solution + solution_actions.base$`1` + solution_actions.base$`2` + solution_actions.base$`3` + solution_actions.base$`4`
#Plot som of actions with tmap library
plot.base <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sum_actions.base", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot.base
The results of this base model will be compared with the subsequent models generated.
This model differs from the previous one in that it tries to group conservation actions within the selected sites as part of the management plan (through a non-linear relationship in the calculation of benefits). The latter is done by adding a value other than 1 to the curve parameter. Where: (1) indicates that there is a linear relationship between the quotient between the actions carried out with respect to the possible actions to be carried out with respect to the benefit obtained by a characteristic. (2) indicates a quadratic relationship between these values, and (3) a cubic relationship. In simple terms, the cubic relationship further penalizes not performing all conservation actions at a site for a particular feature.
input_data.curve <- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_binary_large,
threats = threats_data_binary_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Correctly loaded inputs
model.curve <- prioriactions::min_costs(input_data.curve, blm = 0, blm_actions = 0, curve = 3)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.curve <- prioriactions::solve(model.curve, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 301479 rows, 360551 columns and 1153705 nonzeros
#> Model fingerprint: 0x4d5d908f
#> Variable types: 199412 continuous, 161139 integer (161139 binary)
#> Coefficient statistics:
#> Matrix range [4e-02, 4e+00]
#> Objective range [1e+00, 1e+00]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+00, 2e+02]
#> Presolve removed 49853 rows and 49888 columns (presolve time = 1s) ...
#> Presolve removed 68004 rows and 68039 columns (presolve time = 2s) ...
#> Presolve removed 68004 rows and 74124 columns (presolve time = 3s) ...
#> Presolve removed 68004 rows and 118792 columns (presolve time = 6s) ...
#> Presolve removed 68004 rows and 152707 columns (presolve time = 7s) ...
#> Presolve removed 140554 rows and 152707 columns (presolve time = 7s) ...
#> Presolve removed 140554 rows and 152707 columns (presolve time = 8s) ...
#> Presolve removed 140554 rows and 152707 columns
#> Presolve time: 8.18s
#> Presolved: 160925 rows, 207844 columns, 741013 nonzeros
#> Variable types: 120829 continuous, 87015 integer (87015 binary)
#>
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#>
#>
#> Root simplex log...
#>
#> Iteration Objective Primal Inf. Dual Inf. Time
#> 0 8.7000000e+01 5.665120e+03 8.115843e+10 10s
#> 10566 1.6493948e+03 8.977190e+02 3.622366e+10 11s
#> 15871 1.6847355e+03 6.464839e+02 1.832565e+10 12s
#> 21047 1.7383311e+03 3.691404e+02 1.022477e+10 13s
#> 23073 1.7489576e+03 2.919893e+02 1.960043e+10 14s
#> 24993 1.7558028e+03 2.394743e+02 1.106023e+10 15s
#> 27259 1.7580825e+03 1.751924e+02 4.865399e+09 16s
#> 30577 1.7573666e+03 1.397468e+02 4.231776e+10 17s
#> 33064 1.7615242e+03 1.004403e+02 8.089102e+09 18s
#> 35510 1.7729093e+03 6.790741e+01 4.128746e+08 19s
#> 38664 1.7298006e+03 5.834824e+01 5.082061e+10 20s
#> 40759 1.7414736e+03 4.138714e+01 1.519211e+10 21s
#> 43028 1.7471199e+03 1.674000e+01 1.187171e+10 22s
#> 44854 1.7432058e+03 7.907891e+00 2.539410e+09 23s
#> 46708 1.7129775e+03 1.846839e-03 1.070016e+07 24s
#> 46899 4.0241344e+03 0.000000e+00 1.237047e+04 24s
#> 51428 2.6521439e+03 0.000000e+00 3.729628e+04 25s
#> 54201 2.4894892e+03 0.000000e+00 3.926573e+05 26s
#> 57678 2.4029214e+03 0.000000e+00 1.691656e+05 27s
#> 59896 2.3584194e+03 0.000000e+00 1.116342e+04 28s
#> 62449 2.2881074e+03 0.000000e+00 3.687678e+04 29s
#> 64590 2.2361691e+03 0.000000e+00 2.527561e+04 30s
#> 67292 2.1767205e+03 0.000000e+00 5.255631e+04 31s
#> 69685 2.1264155e+03 0.000000e+00 2.306073e+04 32s
#> 72755 2.0664684e+03 0.000000e+00 2.686298e+04 33s
#> 74663 2.0213556e+03 0.000000e+00 6.719329e+04 34s
#> 77022 2.0004158e+03 0.000000e+00 3.404845e+04 35s
#> 79202 1.9625694e+03 0.000000e+00 1.404740e+05 36s
#> 81281 1.9404075e+03 0.000000e+00 6.751236e+04 37s
#> 83370 1.9075073e+03 0.000000e+00 2.233910e+05 38s
#> 85196 1.8781112e+03 0.000000e+00 3.072513e+04 39s
#> 86778 1.8623475e+03 0.000000e+00 5.848692e+04 40s
#> 88779 1.8478615e+03 0.000000e+00 1.581817e+04 41s
#> 91034 1.8202979e+03 0.000000e+00 4.674761e+04 42s
#> 92316 1.8039199e+03 0.000000e+00 7.399480e+04 43s
#> 93842 1.7893068e+03 0.000000e+00 4.079724e+04 44s
#> 95507 1.7741819e+03 0.000000e+00 9.531235e+04 45s
#> 97406 1.7601923e+03 0.000000e+00 1.973123e+05 46s
#> 98853 1.7457871e+03 0.000000e+00 6.493497e+04 47s
#> 100751 1.7274354e+03 0.000000e+00 6.061209e+04 48s
#> 102767 1.7068078e+03 0.000000e+00 8.597032e+04 49s
#> 104235 1.6937678e+03 0.000000e+00 3.242944e+04 50s
#> 105359 1.6845971e+03 0.000000e+00 1.806713e+04 51s
#> 106915 1.6700077e+03 0.000000e+00 1.898088e+04 52s
#> 109377 1.6536675e+03 0.000000e+00 1.706546e+04 53s
#> 110817 1.6459181e+03 0.000000e+00 3.270242e+04 54s
#> 112683 1.6328625e+03 0.000000e+00 9.579376e+04 55s
#> 113847 1.6228827e+03 0.000000e+00 1.734708e+04 56s
#> 115749 1.6072338e+03 0.000000e+00 8.961769e+04 57s
#> 117172 1.5951116e+03 0.000000e+00 2.660406e+04 58s
#> 118800 1.5815913e+03 0.000000e+00 2.249393e+04 59s
#> 120561 1.5666350e+03 0.000000e+00 2.340476e+04 60s
#> 122398 1.5540525e+03 0.000000e+00 5.955468e+04 61s
#> 123646 1.5478251e+03 0.000000e+00 3.784425e+04 62s
#> 124878 1.5411165e+03 0.000000e+00 7.520535e+04 63s
#> 126638 1.5301915e+03 0.000000e+00 3.504995e+04 64s
#> 128054 1.5241642e+03 0.000000e+00 2.780520e+05 65s
#> 129014 1.5205610e+03 0.000000e+00 5.458844e+04 66s
#> Concurrent spin time: 0.00s
#>
#> Solved with dual simplex
#>
#> Root relaxation: objective 1.039718e+03, 30576 iterations, 57.01 seconds
#> Total elapsed time = 70.58s
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 1039.71757 0 11890 - 1039.71757 - - 73s
#> 0 0 1186.30479 0 8871 - 1186.30479 - - 139s
#> 0 0 1192.02160 0 8333 - 1192.02160 - - 174s
#> 0 0 1192.62345 0 8579 - 1192.62345 - - 179s
#> 0 0 1192.72178 0 8479 - 1192.72178 - - 182s
#> 0 0 1192.73665 0 8435 - 1192.73665 - - 184s
#> 0 0 1252.04391 0 6723 - 1252.04391 - - 257s
#> 0 0 1254.79649 0 6584 - 1254.79649 - - 266s
#> 0 0 1254.97535 0 6618 - 1254.97535 - - 270s
#> 0 0 1255.01085 0 6658 - 1255.01085 - - 275s
#> 0 0 1304.65302 0 5183 - 1304.65302 - - 414s
#> 0 0 1319.82558 0 5053 - 1319.82558 - - 449s
#> 0 0 1320.37935 0 5093 - 1320.37935 - - 465s
#> 0 0 1320.49364 0 5086 - 1320.49364 - - 476s
#> 0 0 1320.53796 0 5055 - 1320.53796 - - 491s
#> 0 0 1362.32645 0 3870 - 1362.32645 - - 588s
#> H 0 0 1896.0000000 1362.32645 28.1% - 589s
#> H 0 0 1860.0000000 1362.32645 26.8% - 589s
#> 0 0 1379.11885 0 3223 1860.00000 1379.11885 25.9% - 651s
#> 0 0 1380.24570 0 3094 1860.00000 1380.24570 25.8% - 674s
#> 0 0 1380.31869 0 3156 1860.00000 1380.31869 25.8% - 688s
#> 0 0 1405.41632 0 2114 1860.00000 1405.41632 24.4% - 771s
#> H 0 0 1624.0000000 1405.41632 13.5% - 774s
#>
#> Cutting planes:
#> Gomory: 191
#> Cover: 38012
#> Implied bound: 59
#> Clique: 199
#> MIR: 5357
#> Flow cover: 3408
#> RLT: 3907
#> Relax-and-lift: 59
#>
#> Explored 1 nodes (128733 simplex iterations) in 774.63 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 3: 1624 1860 1896
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.624000000000e+03, best bound 1.406000000000e+03, gap 13.4236%
Note that the new model increased its dimensions compared to the previous one (from 11580 variables and 2361 constraints to 207844 variables and 160925 constraints), and therefore, increasing its resolution complexity. In this way, we obtain a value of the objective function of 1497 (greater than the previous one) with a quality of 4.48%.
#Getting unit distribution selected
solution_units.curve <- solution.curve$getSolutionUnits()
#Getting action distribution
solution_actions.curve <- solution.curve$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.curve <- solution_actions.curve$`1`
shp_mitchell$action_2.curve <- solution_actions.curve$`2`
shp_mitchell$action_3.curve <- solution_actions.curve$`3`
shp_mitchell$action_4.curve <- solution_actions.curve$`4`
shp_mitchell$sum_actions.curve <- solution_units.curve$solution + solution_actions.curve$`1` + solution_actions.curve$`2` + solution_actions.curve$`3` + solution_actions.curve$`4`
#Plot som of actions with tmap library
plot.curve <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sum_actions.curve", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) +
tmap::tm_borders(col="black", lwd = 0.5)
#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.curve)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Thus, we can notice that there are a greater number of sites with a higher density of actions with respect to the previous model (for example, a greater number of units painted yellow).
To add spatial requirements to the model, there are two key parameters: blm which works in the same way as the parameter of the same name in marxan, which tries to minimize fragmentation between the selected planning units (i.e. regardless if conservation actions are carried out within them). Next we will see what happens when this parameter is relevant in the creation of the model.
input_data.spatial_req<- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_binary_large,
threats = threats_data_binary_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Correctly loaded inputs
model.spatial_req <- prioriactions::min_costs(input_data.spatial_req, blm = 10, blm_actions = 0, curve = 1)
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.spatial_req <- prioriactions::solve(model.spatial_req, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 43557 rows, 25312 columns and 202916 nonzeros
#> Model fingerprint: 0x3fe989f5
#> Variable types: 0 continuous, 25312 integer (25312 binary)
#> Coefficient statistics:
#> Matrix range [3e-01, 4e+00]
#> Objective range [1e+00, 1e+02]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+00, 2e+02]
#> Found heuristic solution: objective 3304.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.21s
#> Presolved: 43557 rows, 25277 columns, 202881 nonzeros
#> Variable types: 0 continuous, 25277 integer (25277 binary)
#>
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#>
#>
#> Root simplex log...
#>
#> Iteration Objective Primal Inf. Dual Inf. Time
#> 31094 3.0518120e+03 0.000000e+00 1.657631e+05 1s
#> 36117 2.2436091e+03 0.000000e+00 2.109409e+05 2s
#> 40599 2.1094886e+03 0.000000e+00 1.042250e+05 3s
#> 44622 1.9300104e+03 0.000000e+00 4.997226e+04 4s
#> 46823 1.8125023e+03 0.000000e+00 5.160657e+06 5s
#> 48643 1.7211745e+03 0.000000e+00 7.747742e+04 6s
#> 50285 1.6811883e+03 0.000000e+00 1.170770e+07 7s
#> 51800 1.6155709e+03 0.000000e+00 2.597011e+05 8s
#> 53093 1.5878591e+03 0.000000e+00 5.304713e+04 9s
#> Concurrent spin time: 0.60s
#>
#> Solved with dual simplex
#>
#> Root relaxation: objective 1.094585e+03, 15227 iterations, 9.03 seconds
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 1094.58473 0 12761 3304.00000 1094.58473 66.9% - 10s
#> H 0 0 3296.0000000 1094.58473 66.8% - 11s
#> 0 0 1243.82454 0 18610 3296.00000 1243.82454 62.3% - 29s
#> 0 0 1490.21388 0 13269 3296.00000 1490.21388 54.8% - 43s
#> 0 0 1490.21469 0 13271 3296.00000 1490.21469 54.8% - 56s
#> 0 0 1490.21469 0 13272 3296.00000 1490.21469 54.8% - 70s
#> 0 0 1490.21469 0 13273 3296.00000 1490.21469 54.8% - 70s
#> 0 0 1490.21469 0 13269 3296.00000 1490.21469 54.8% - 80s
#> H 0 0 2414.0927254 1490.21469 38.3% - 83s
#> H 0 0 2267.1903810 1490.21469 34.3% - 84s
#> H 0 0 2245.1903810 1490.21469 33.6% - 84s
#> H 0 2 2235.1903810 1490.21469 33.3% - 85s
#> 0 2 1490.21469 0 13267 2235.19038 1490.21469 33.3% - 85s
#> 1 4 1506.90672 1 13296 2235.19038 1490.24408 33.3% 5636 99s
#> 3 6 1520.70894 2 12947 2235.19038 1507.43883 32.6% 7508 115s
#> 5 8 1596.13436 2 13427 2235.19038 1520.76378 32.0% 7228 127s
#> 7 10 1605.31844 3 12311 2235.19038 1520.78909 32.0% 6821 134s
#> 9 12 1530.71185 3 12332 2235.19038 1530.71185 31.5% 5688 211s
#> 11 14 1551.48967 4 12282 2235.19038 1530.71351 31.5% 5743 222s
#> 13 16 1575.02183 5 12700 2235.19038 1530.71351 31.5% 5761 230s
#> 15 18 1622.37087 5 12084 2235.19038 1530.71351 31.5% 5646 238s
#> 17 20 1628.58446 6 11983 2235.19038 1530.71351 31.5% 5651 249s
#> 19 23 1579.58649 6 12678 2235.19038 1530.71351 31.5% 5598 263s
#> 22 25 1620.57003 8 12483 2235.19038 1530.71351 31.5% 5930 275s
#> 24 27 1595.54701 8 12536 2235.19038 1530.71351 31.5% 5901 286s
#> H 26 29 2233.1903810 1530.71351 31.5% 5689 301s
#> H 27 29 2225.1903810 1530.71351 31.2% 5872 301s
#> 28 31 1644.45520 9 11236 2225.19038 1530.71351 31.2% 6237 316s
#> H 29 31 2193.8795220 1530.71351 30.2% 6021 316s
#> 30 33 1626.77311 10 12249 2193.87952 1530.71351 30.2% 6292 327s
#> 32 35 1655.34224 11 12501 2193.87952 1530.71351 30.2% 6254 337s
#> 34 37 1674.91212 11 10028 2193.87952 1530.71351 30.2% 6394 348s
#> 36 40 1701.99427 12 8894 2193.87952 1530.71351 30.2% 6414 359s
#> 39 44 1735.29735 13 10116 2193.87952 1530.71351 30.2% 6320 375s
#> 43 51 1726.77391 14 7611 2193.87952 1530.71351 30.2% 6189 388s
#> 50 55 1764.48171 18 6790 2193.87952 1530.71351 30.2% 5862 400s
#> H 54 61 2172.8795220 1530.71351 29.6% 5816 412s
#> 60 69 1779.97493 21 5736 2172.87952 1530.71351 29.6% 5655 425s
#> 68 74 1809.64250 25 3476 2172.87952 1530.71351 29.6% 5321 436s
#> 73 87 1821.36263 26 2316 2172.87952 1530.71351 29.6% 5234 447s
#> H 86 99 1845.0253024 1530.71351 17.0% 4652 454s
#> H 90 99 1824.1807319 1530.71351 16.1% 4447 454s
#> H 92 99 1823.1807319 1530.71351 16.0% 4351 455s
#>
#> Cutting planes:
#> Zero half: 1
#>
#> Explored 98 nodes (452968 simplex iterations) in 455.03 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 10: 1823.18 1824.18 1845.03 ... 2267.19
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.823180731900e+03, best bound 1.530713509025e+03, gap 16.0416%
In the same way as what happened in the previous model, it is naturally inferred that the spatial requirement implies higher costs in the management plan, as well as a model with greater computational complexity to resolve.
#Getting unit distribution selected
solution_units.spatial_req <- solution.spatial_req$getSolutionUnits()
#Getting action distribution
solution_actions.spatial_req<- solution.spatial_req$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.spatial_req <- solution_actions.spatial_req$`1`
shp_mitchell$action_2.spatial_req <- solution_actions.spatial_req$`2`
shp_mitchell$action_3.spatial_req <- solution_actions.spatial_req$`3`
shp_mitchell$action_4.spatial_req <- solution_actions.spatial_req$`4`
shp_mitchell$sum_actions.spatial_req <- solution_units.spatial_req$solution + solution_actions.spatial_req$`1` + solution_actions.spatial_req$`2` + solution_actions.spatial_req$`3` + solution_actions.spatial_req$`4`
#Plot som of actions with tmap library
plot.spatial_req <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sum_actions.spatial_req", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) +
tmap::tm_borders(col="black", lwd = 0.5)
#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.spatial_req)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
As seen, minimizing fragmentation produced a more closely related conservation plan compared to the base model.
Another parameter to take into account is the blm_actions that adds the minimization of fragmentation between conservation actions, this means that it tries to spatially bring the management actions closer together and therefore, independently achieve a management plan as well related.
input_data.spatial_req_actions<- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_binary_large,
threats = threats_data_binary_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Correctly loaded inputs
model.spatial_req_actions <- prioriactions::min_costs(input_data.spatial_req_actions, blm = 0, blm_actions = 10, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
solution.spatial_req_actions <- prioriactions::solve(model.spatial_req_actions, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 167145 rows, 66508 columns and 491288 nonzeros
#> Model fingerprint: 0x8086a2c3
#> Variable types: 0 continuous, 66508 integer (66508 binary)
#> Coefficient statistics:
#> Matrix range [3e-01, 4e+00]
#> Objective range [1e+00, 1e+02]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+00, 2e+02]
#> Found heuristic solution: objective 11548.000000
#> Presolve time: 0.83s
#> Presolved: 167145 rows, 66508 columns, 491288 nonzeros
#> Variable types: 0 continuous, 66508 integer (66508 binary)
#>
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#>
#>
#> Root simplex log...
#>
#> Iteration Objective Primal Inf. Dual Inf. Time
#> 0 4.4417447e+04 0.000000e+00 4.631422e+05 2s
#> Warning: Markowitz tolerance tightened to 0.5
#> 61617 2.4718346e+04 1.030990e-03 2.720499e+07 2s
#> 66858 2.2891668e+04 0.000000e+00 3.852743e+05 3s
#> 69133 2.1911720e+04 0.000000e+00 2.224889e+06 4s
#> 71362 2.1302148e+04 0.000000e+00 4.572488e+05 5s
#> 73592 2.0670783e+04 0.000000e+00 1.453689e+06 6s
#> 75648 2.0151453e+04 0.000000e+00 4.117813e+05 7s
#> 77622 1.9833790e+04 0.000000e+00 4.677604e+05 8s
#> 79253 1.9750604e+04 0.000000e+00 4.578073e+05 9s
#> 80145 1.8872755e+04 0.000000e+00 1.105365e+07 10s
#> 81929 1.7169011e+04 0.000000e+00 1.248981e+07 11s
#> 83267 1.6652827e+04 0.000000e+00 9.086708e+05 12s
#> 84605 1.5942891e+04 0.000000e+00 1.345409e+06 13s
#> 86166 1.5443126e+04 0.000000e+00 5.417270e+06 14s
#> 87504 1.5222457e+04 0.000000e+00 2.737745e+07 15s
#> 88620 1.4570497e+04 0.000000e+00 1.905453e+06 16s
#> 89958 1.4107801e+04 0.000000e+00 1.628268e+06 17s
#> 91296 1.3683932e+04 0.000000e+00 2.214303e+06 18s
#> 92857 1.3290566e+04 0.000000e+00 2.393024e+06 19s
#> 94195 1.3000940e+04 0.000000e+00 4.811312e+06 20s
#> 95533 1.2857727e+04 0.000000e+00 1.515921e+06 21s
#> 97094 1.2696943e+04 0.000000e+00 1.366198e+06 22s
#> 98655 1.2491999e+04 0.000000e+00 2.139967e+06 23s
#> 99547 1.2379188e+04 0.000000e+00 2.647626e+06 24s
#> 100885 1.2050869e+04 0.000000e+00 5.010568e+06 25s
#> 102446 1.1910771e+04 0.000000e+00 1.677002e+06 26s
#> 104007 1.1800496e+04 0.000000e+00 1.166538e+06 27s
#> 105791 1.1582139e+04 0.000000e+00 1.531039e+06 28s
#> 107129 1.1310217e+04 0.000000e+00 4.454977e+05 29s
#> 108467 1.1105586e+04 0.000000e+00 9.934501e+05 30s
#> 110028 1.0944070e+04 0.000000e+00 4.081888e+05 31s
#> 111143 1.0860066e+04 0.000000e+00 8.532527e+07 32s
#> 112704 1.0211889e+04 0.000000e+00 5.370109e+06 33s
#> 114042 9.8247835e+03 0.000000e+00 2.867375e+06 34s
#> Concurrent spin time: 0.00s
#>
#> Solved with dual simplex
#>
#> Root relaxation: objective 1.289311e+03, 19774 iterations, 33.49 seconds
#> Total elapsed time = 35.10s
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 1289.31051 0 17015 11548.0000 1289.31051 88.8% - 40s
#> H 0 0 10906.955285 1289.31051 88.2% - 41s
#> H 0 0 9770.6621484 1289.31051 86.8% - 41s
#> 0 0 1322.21108 0 28897 9770.66215 1322.21108 86.5% - 76s
#> 0 0 1691.79293 0 31021 9770.66215 1691.79293 82.7% - 169s
#> H 0 0 5843.0029202 1691.79293 71.0% - 169s
#> H 0 0 5842.0029202 1691.79293 71.0% - 169s
#> H 0 0 5736.9354972 1691.79293 70.5% - 169s
#> 0 0 1691.79293 0 36204 5736.93550 1691.79293 70.5% - 225s
#> H 0 0 5486.5993358 1691.79293 69.2% - 225s
#> H 0 0 5464.3631894 1691.79293 69.0% - 226s
#> 0 0 1691.82283 0 36205 5464.36319 1691.82283 69.0% - 226s
#> 0 0 1691.82283 0 36205 5464.36319 1691.82283 69.0% - 304s
#> H 0 0 5309.6680442 1691.87933 68.1% - 338s
#> H 0 0 4213.3417431 1691.87933 59.8% - 342s
#> H 0 2 4213.1863136 1691.87933 59.8% - 343s
#> 0 2 1691.87933 0 36205 4213.18631 1691.87933 59.8% - 343s
#> 1 4 1709.96480 1 36294 4213.18631 1691.87933 59.8% 3309 365s
#> 3 6 1730.99919 2 36432 4213.18631 1709.96695 59.4% 3502 393s
#> 5 8 1736.65762 2 26842 4213.18631 1721.29090 59.1% 3103 409s
#> 7 10 1764.96275 3 30424 4213.18631 1731.00579 58.9% 3013 421s
#> 9 12 1758.68503 3 26930 4213.18631 1731.00579 58.9% 2860 435s
#> 11 14 1771.40852 4 23138 4213.18631 1731.00579 58.9% 2831 454s
#> 13 16 1793.83908 5 23135 4213.18631 1731.00579 58.9% 3002 465s
#> 15 18 1832.32558 5 22295 4213.18631 1731.00579 58.9% 2861 476s
#> 17 20 1841.19368 6 21344 4213.18631 1731.00579 58.9% 2845 492s
#> 19 22 1810.76823 6 25890 4213.18631 1731.00579 58.9% 2804 504s
#> 21 24 1846.91216 7 22504 4213.18631 1731.00579 58.9% 2701 515s
#> 23 26 1860.50029 7 23807 4213.18631 1731.00579 58.9% 2653 525s
#> 25 29 1847.60307 8 22424 4213.18631 1731.00579 58.9% 2555 532s
#> H 28 31 3496.0563978 1731.00579 50.5% 2363 555s
#> H 29 31 3199.6890502 1731.00579 45.9% 2324 555s
#> 30 35 1868.08674 10 22671 3199.68905 1731.00579 45.9% 2461 570s
#> 34 37 1885.18418 10 20161 3199.68905 1731.00579 45.9% 2350 586s
#> 36 41 1894.64458 11 21364 3199.68905 1731.00579 45.9% 2352 606s
#> 40 45 1894.60195 12 19349 3199.68905 1731.00579 45.9% 2306 627s
#> 44 49 1920.74844 14 20776 3199.68905 1731.00579 45.9% 2285 645s
#> 48 55 1927.96817 15 20650 3199.68905 1731.00579 45.9% 2243 660s
#> 54 59 1971.17510 18 20518 3199.68905 1731.00579 45.9% 2144 670s
#> H 58 63 3097.7023441 1731.00579 44.1% 2045 686s
#> H 61 63 2538.3060493 1731.00579 31.8% 2002 686s
#> 62 67 2004.85523 21 20189 2538.30605 1731.00579 31.8% 2024 698s
#> 66 72 2021.82057 22 17772 2538.30605 1731.00579 31.8% 1996 711s
#> 71 75 2033.33996 23 20909 2538.30605 1731.00579 31.8% 1980 768s
#> 74 81 2073.08900 24 23025 2538.30605 1731.00579 31.8% 1959 783s
#> 80 88 2005.55142 25 20002 2538.30605 1731.00579 31.8% 1892 797s
#> 87 92 2014.19431 26 20619 2538.30605 1731.00579 31.8% 1826 822s
#> H 88 92 2301.6723228 1731.00579 24.8% 1805 822s
#> 91 101 2014.51546 27 19985 2301.67232 1731.00579 24.8% 1783 837s
#> 100 110 2069.48068 30 23060 2301.67232 1731.00579 24.8% 1756 858s
#> 109 119 2064.17616 31 23141 2301.67232 1731.00579 24.8% 1724 875s
#> 118 133 2070.25107 36 17861 2301.67232 1731.00579 24.8% 1668 891s
#> 132 144 2098.22708 41 15863 2301.67232 1731.00579 24.8% 1581 915s
#> 143 163 2122.03727 45 12823 2301.67232 1731.00579 24.8% 1552 942s
#> H 154 163 2277.9698879 1731.00579 24.0% 1499 942s
#> 162 176 2165.89784 50 12909 2277.96989 1731.00579 24.0% 1471 961s
#> 175 186 2201.05396 54 15225 2277.96989 1731.00579 24.0% 1453 986s
#> 189 188 2247.34693 59 17095 2277.96989 1731.00579 24.0% 1437 1013s
#> 203 191 2268.44703 63 13188 2277.96989 1746.19068 23.3% 1387 1034s
#> 214 196 1756.58558 4 35918 2277.96989 1746.19677 23.3% 1358 1061s
#> 219 201 1786.77530 6 35835 2277.96989 1746.19677 23.3% 1374 1097s
#> 224 205 1812.63049 7 33333 2277.96989 1746.19677 23.3% 1394 1125s
#> 228 209 1835.56037 9 32971 2277.96989 1746.19677 23.3% 1413 1159s
#> H 232 202 2261.9698879 1746.19677 22.8% 1425 1194s
#> 236 209 1875.26386 12 30344 2261.96989 1746.19677 22.8% 1458 1226s
#> 243 216 1938.96010 13 23585 2261.96989 1746.19677 22.8% 1484 1253s
#> 250 228 1959.08680 15 23496 2261.96989 1746.19677 22.8% 1484 1281s
#> H 262 238 2258.7204477 1746.19677 22.7% 1467 1322s
#> 273 251 1988.74581 18 27565 2258.72045 1746.19677 22.7% 1473 1352s
#> 286 267 2012.31352 21 29597 2258.72045 1746.19677 22.7% 1465 1386s
#> 302 281 2075.94227 22 12576 2258.72045 1746.19677 22.7% 1460 1421s
#> 316 310 2094.87900 27 15761 2258.72045 1746.19677 22.7% 1472 1455s
#> 345 344 2144.34428 28 20344 2258.72045 1746.19677 22.7% 1428 1499s
#> H 349 334 2242.1120235 1746.19677 22.1% 1421 1499s
#> 381 333 2132.12647 34 19645 2242.11202 1746.19677 22.1% 1347 1537s
#> 394 339 2149.27546 39 22322 2242.11202 1746.19677 22.1% 1339 1634s
#> H 397 254 2165.5783120 1746.19677 19.4% 1331 1634s
#>
#> Cutting planes:
#> Zero half: 1
#>
#> Explored 401 nodes (594133 simplex iterations) in 1634.78 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 10: 2165.58 2242.11 2258.72 ... 3496.06
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 2.165578312040e+03, best bound 1.746196772342e+03, gap 19.3658%
The effect of incorporating this parameter (blm_actions) can be seen better in the action distributions:
#Getting unit distribution selected
solution_units.spatial_req_actions <- solution.spatial_req_actions$getSolutionUnits()
#Getting action distribution
solution_actions.spatial_req_actions<- solution.spatial_req_actions$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.spatial_req_actions <- solution_actions.spatial_req_actions$`1`
shp_mitchell$action_2.spatial_req_actions <- solution_actions.spatial_req_actions$`2`
shp_mitchell$action_3.spatial_req_actions <- solution_actions.spatial_req_actions$`3`
shp_mitchell$action_4.spatial_req_actions <- solution_actions.spatial_req_actions$`4`
#Actions plots
plot_action1.spatial_req_actions <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_1.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action2.spatial_req_actions <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_2.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action3.spatial_req_actions <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_3.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
plot_action4.spatial_req_actions <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("action_4.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) +
tmap::tm_borders(col="black", lwd = 0.5)
tmap::tmap_arrange(plot_action1.spatial_req_actions , plot_action2.spatial_req_actions , plot_action3.spatial_req_actions , plot_action4.spatial_req_actions)
An interesting comparison is how the use of continuous threat intensities influences solutions with respect to the use of binary threat presence / absence data. To do this, we change the column amount for the original output from the simulation of threat distribution, without performing a previous filter, i.e.,
#We get the extended matrix of distribution data for both cases, set the specific col names and we are left with only positive values
#threat distribution data-------------------------------------
threats_data_large <- reshape2::melt(threats_data)
colnames(threats_data_large) <- c("pu", "threats", "amount")
threats_data_large <- threats_data_large[threats_data_large$amount > 0, ]
threats_data_large$cost <- 1
threats_data_large$status <- 0
Thus, the new model is build as follow:
input_data.threat_cont <- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_binary_large,
threats = threats_data_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Correctly loaded inputs
model.threat_cont<- prioriactions::min_costs(input_data.threat_cont, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.threat_cont <- prioriactions::solve(model.threat_cont, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 106792 nonzeros
#> Model fingerprint: 0xf3148484
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#> Matrix range [3e-02, 4e+00]
#> Objective range [1e+00, 1e+00]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+02, 2e+02]
#> Found heuristic solution: objective 3203.0000000
#> Found heuristic solution: objective 1934.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.08s
#> Presolved: 2361 rows, 11545 columns, 106757 nonzeros
#> Variable types: 0 continuous, 11545 integer (11545 binary)
#>
#> Root relaxation: objective 7.061076e+02, 4161 iterations, 0.28 seconds
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 706.10764 0 596 1934.00000 706.10764 63.5% - 0s
#> H 0 0 1142.0000000 706.10764 38.2% - 0s
#> H 0 0 1131.0000000 706.10764 37.6% - 0s
#> 0 0 832.81049 0 742 1131.00000 832.81049 26.4% - 0s
#> H 0 0 1102.0000000 832.81049 24.4% - 1s
#> 0 0 832.84104 0 742 1102.00000 832.84104 24.4% - 1s
#> 0 0 932.02978 0 621 1102.00000 932.02978 15.4% - 1s
#>
#> Cutting planes:
#> Cover: 1152
#> Clique: 4
#> MIR: 24
#>
#> Explored 1 nodes (7980 simplex iterations) in 1.58 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 5: 1102 1131 1142 ... 3203
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.102000000000e+03, best bound 9.330000000000e+02, gap 15.3358%
#Getting unit distribution selected
solution_units.threat_cont <- solution.threat_cont$getSolutionUnits()
#Getting action distribution
solution_actions.threat_cont <- solution.threat_cont$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.threat_cont <- solution_actions.threat_cont$`1`
shp_mitchell$action_2.threat_cont <- solution_actions.threat_cont$`2`
shp_mitchell$action_3.threat_cont <- solution_actions.threat_cont$`3`
shp_mitchell$action_4.threat_cont <- solution_actions.threat_cont$`4`
shp_mitchell$sum_actions.threat_cont <- solution_units.threat_cont$solution + solution_actions.threat_cont$`1` + solution_actions.threat_cont$`2` + solution_actions.threat_cont$`3` + solution_actions.threat_cont$`4`
#Plot som of actions with tmap library
plot.threat_cont <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sum_actions.threat_cont", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) +
tmap::tm_borders(col="black", lwd = 0.5)
#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.threat_cont)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
The more general model allows us to work with continuous data on intensity of species and threats. For this, the same procedure is carried out on the data done in the previous step, but now with the amount column of the features. Also, it is necessary to update the target column from the features_data file.
#Feature distribution data-------------------------------------
rij_data_large <- reshape2::melt(rij_data)
colnames(rij_data_large) <- c("pu", "species", "amount")
rij_data_large <- rij_data_large[rij_data_large$amount > 0, ]
#features data-------------------------------------
features_data$target <- colSums(rij_data)*0.15
Thus, the new model is build as follow:
input_data.feature_cont <- prioriactions::problem(
pu = pu_data, features = features_data, rij = rij_data_large,
threats = threats_data_large, sensitivity = sensitivity_data,
bound = bound_data
)
#> Correctly loaded inputs
model.feature_cont<- prioriactions::min_costs(input_data.feature_cont, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.feature_cont <- prioriactions::solve(model.feature_cont, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 213072 nonzeros
#> Model fingerprint: 0x73cee9ae
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#> Matrix range [2e-03, 4e+00]
#> Objective range [1e+00, 1e+00]
#> Bounds range [1e+00, 1e+00]
#> RHS range [1e+02, 2e+02]
#> Found heuristic solution: objective 3302.0000000
#> Found heuristic solution: objective 2299.0000000
#> Presolve time: 0.20s
#> Presolved: 2361 rows, 11580 columns, 213072 nonzeros
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#>
#> Root relaxation: objective 8.270739e+02, 2992 iterations, 0.17 seconds
#>
#> Nodes | Current Node | Objective Bounds | Work
#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
#>
#> 0 0 827.07394 0 662 2299.00000 827.07394 64.0% - 0s
#> H 0 0 1317.0000000 827.07394 37.2% - 0s
#> 0 0 959.14938 0 789 1317.00000 959.14938 27.2% - 1s
#> H 0 0 1289.0000000 959.14938 25.6% - 1s
#> 0 0 959.25344 0 790 1289.00000 959.25344 25.6% - 1s
#> 0 0 959.25585 0 790 1289.00000 959.25585 25.6% - 1s
#> 0 0 1054.37015 0 788 1289.00000 1054.37015 18.2% - 2s
#>
#> Cutting planes:
#> Cover: 1323
#> Clique: 2
#> MIR: 39
#>
#> Explored 1 nodes (6483 simplex iterations) in 2.29 seconds
#> Thread count was 2 (of 4 available processors)
#>
#> Solution count 4: 1289 1317 2299 3302
#>
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.289000000000e+03, best bound 1.055000000000e+03, gap 18.1536%
#Getting unit distribution selected
solution_units.feature_cont <- solution.feature_cont$getSolutionUnits()
#Getting action distribution
solution_actions.feature_cont <- solution.feature_cont$getSolutionActions()
#Assign solution to shapefile field to plot it
shp_mitchell$action_1.feature_cont <- solution_actions.feature_cont$`1`
shp_mitchell$action_2.feature_cont <- solution_actions.feature_cont$`2`
shp_mitchell$action_3.feature_cont <- solution_actions.feature_cont$`3`
shp_mitchell$action_4.feature_cont <- solution_actions.feature_cont$`4`
shp_mitchell$sum_actions.feature_cont <- solution_units.feature_cont$solution + solution_actions.feature_cont$`1` + solution_actions.feature_cont$`2` + solution_actions.feature_cont$`3` + solution_actions.feature_cont$`4`
#Plot sum of actions with tmap library
plot.feature_cont <- tmap::tm_shape(shp_mitchell) +
tmap::tm_fill("sum_actions.feature_cont", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) +
tmap::tm_borders(col="black", lwd = 0.5)
#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.feature_cont)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.